Native Plants Analysis

library(vegan)
library(broom)
library(MASS)
library(tidyverse)

## Do Not Run
#wExp <- readRDS("~/RMB/SoilDomestication/Data/wExp.rds")
#wExp %>% 
#  filter(Site == "Arkansas") %>% 
#  ungroup() %>% 
#  saveRDS("~/RMB/SoilDomestication/Data/native_plants_data.rds")
nExp <- readRDS("~/RMB/SoilDomestication/Data/native_plants_data.rds") %>% filter(SampleID != "Arw.22")
tax <- readRDS("~/RMB/SoilDomestication/Data/gg_otus_tax.rds")

long_pcoa <- function(x, samples = "SampleID", otus = "variable", value = "RA", dist = "bray"){
  to_drop <- otus
  metadata <- x %>% 
    dplyr::ungroup() %>% 
    purrr::discard(is.double) %>% 
    dplyr::select_(.dots = paste("-", to_drop)) %>% 
    dplyr::distinct()
  
  wide_table <- x %>% 
    dplyr::select_(samples, otus, value) %>% 
    tidyr::spread_(otus, value, fill = 0)
  
  pc <- vegan::capscale(log2(wide_table[,2:ncol(wide_table)] + 1) ~ 1, dist = dist)
  axes <- dplyr::bind_cols(metadata, dplyr::as_tibble(vegan::scores(pc, choices = c(1:5))$sites))
  eigen_vals <- vegan::eigenvals(pc) / sum(vegan::eigenvals(pc))
  
  return(list(axes = axes, eigen_vals = eigen_vals))
}
raref <- function(x, sampling_depth, value = "value", otus = "variable") {
  rare_values <- data.frame(table(sample(x$`variable`, sampling_depth, replace = T, prob = x$`value`/x$depth)))
  names(rare_values) <- c(otus, "rare_value")
  return(suppressMessages(left_join(x, rare_values, by = otus) %>% replace_na(list(rare_value = 0))))
}
long_adonis <- function(x, samples = "SampleID", otus = "variable", value = "RA", dist = "bray", formula) {
    to_drop <- otus
    metadata <- x %>% 
        dplyr::ungroup() %>% 
        purrr::discard(is.double) %>% 
        dplyr::select_(.dots = paste("-", to_drop)) %>% 
        dplyr::distinct()
  
    wide_table <- x %>% 
        dplyr::select_(samples, otus, value) %>% 
        tidyr::spread_(otus, value, fill = 0)
    permanova <- vegan::adonis(as.formula(paste("wide_table[,2:ncol(wide_table)] ~ ", formula, sep = "")), data = metadata, dist = dist)
    return(permanova)
}

Phyla Stuff

Differential Abundance

pois_glm %>% 
  unnest(map(models, ~tidy(.))) %>% 
  mutate(p.adj = p.adjust(p.value, "fdr")) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(direction = ifelse(estimate > 0, "Weed", "Rice")) %>% 
  filter(p.adj <= 0.05) %>% 
  group_by(Compartment, direction, variable) %>% 
  summarise(n = n()) %>% 
  group_by(Compartment, direction) %>% 
  mutate(total = n()) %>% 
  filter(n == 3) %>% 
  group_by(Compartment, direction, total) %>% 
  summarise(n = n())

nb_glm %>% 
  unnest(map(models, ~tidy(.))) %>% 
  mutate(p.adj = p.adjust(p.value, "fdr")) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(direction = ifelse(estimate > 0, "Weed", "Rice")) %>% 
  filter(p.adj <= 0.05) %>% 
  group_by(Compartment, direction, variable) %>% 
  summarise(n = n()) %>% 
  group_by(Compartment, direction) %>% 
  mutate(total = n()) %>% 
  filter(n == 3) %>% 
  group_by(Compartment, direction, total) %>% 
  summarise(n = n())

plant_lm %>% 
  unnest(map(models, ~tidy(.))) %>% 
  mutate(p.adj = p.adjust(p.value, "fdr")) %>% 
  filter(term != "(Intercept)") %>% 
  mutate(direction = ifelse(estimate > 0, "Weed", "Rice")) %>% 
  filter(p.adj <= 0.05) %>% 
  group_by(Compartment, direction, variable) %>% 
  summarise(n = n()) %>% 
  group_by(Compartment, direction) %>% 
  mutate(total = n()) %>% 
  filter(n == 3) %>% 
  group_by(Compartment, direction, total) %>% 
  summarise(n = n())

nExp2 %>% 
  filter(variable == "1110462" & Compartment == "Rhizosphere") %>% 
  ggplot(aes(host_common_name, rare_value, color = host_common_name)) +
  geom_jitter(height = 0, width = 0.1)
tax["1091625",]